Shear Modulus Estimation by Application of Spatially Modulated Impulse Acoustic Radiation Force Approximation

ABSTRACT

A method for determining a shear modulus of an elastic material with a known density value is provided. In this method, a spatially modulated acoustic radiation force is used to initially generate a disturbance of known spatial frequency or wavelength. The propagation of this initial displacement as a shear wave is measured using ultrasound tracking methods. A temporal frequency is determined based on the shear wave. The shear modulus of the elastic material at the point of excitation may be calculated using the values of the spatial wavelength, material density, and temporal frequency.

REFERENCE TO RELATED APPLICATION

The present application is a continuation of U.S. patent applicationSer. No. 12/118,359, filed May 9, 2008, which in turn claims the benefitof U.S. Provisional Application No. 60/917,031, filed May 9, 2007; thedisclosures of each are hereby incorporated by reference in theirentirety into the present application.

FIELD OF THE INVENTION

The present invention relates to a method for determining a shearmodulus, i.e., stiffness, of an elastic material such as human or animalsoft tissues. More particularly, the method of the present inventionprovides an ultrasonic measurement of tissue shear modulus for thediagnosis of disease, such as liver fibrosis and prostate cancer.

DESCRIPTION OF THE RELATED ART

Quantification of tissue shear modulus (“stiffness”) is needed fordisease diagnosis, therapy monitoring, and biomechanics research.Changes in tissue stiffness have long been associated with disease orinjury. Manual palpation has been used to detect these changes. Manualpalpation is limited to superficial or accessible tissues. Nonetheless,manual palpation can reveal structures that may not be evident inconventional imaging modalities. For instance, prostate cancer is ofteninvisible to ultrasound, while the stiffening of the prostate isdetectable by palpation. Radio frequency (RF) and high intensity focusedultrasound (HIFU) ablation can cause tissue necrosis and stiffening oftissue, while being invisible to ultrasound.

The limitations and diagnostic utility of palpation have led to aninterest in elastography, the generation of images of or related totissue stiffness. These images are useful because disease or therapeuticprocesses may cause a significant change in tissue stiffness without aconcomitant change in ultrasound echogenicity or x-ray density. A numberof techniques have been explored including strain elastography,transient elastography, sonoelastography, vibroacoustograpy, andacoustic radiation force impulse imaging (ARFI). ARFI is attractivebecause of potential to induce measureable displacement wherever anultrasound beam may penetrate. Sonoelasticity methods, both transientand steady-state have the potential to provide a quantitative estimateof shear modulus.

It is well known that the shear modulus G of a linear elastic materialwith density ρ is related to the speed of shear wave propagation c_(s)by c_(s)=√{square root over (G/ρ)}. Many authors have used thisrelationship to estimate the local shear modulus of tissues. In one ofthese methods, a point or extended source is driven at a givenfrequency, and the local wavelength of the shear wave is estimated byultrasonic or magnetic resonance imaging (MRI) measurements of theinduced displacement field. A variety of techniques are employed but alldepend on a source of known temporal frequency. Other methods employstep excitation to generate non-harmonic shear waves, and local shearmodulus properties are estimated by inversion of the wave equation orarrival time methods. However, the inversion method is sensitive tonoise.

The background theory of the present invention is as follows:

In a linear elastic material the shear component of wave motion isgoverned by the homogenous wave equation:

$\begin{matrix}{{{\rho \; \frac{\partial^{2}u_{i}}{\partial t^{2}}} - {G\; {\nabla^{2}u_{i}}}} = 0} & (1)\end{matrix}$

where μ_(i) is the displacement in the i direction, G is the shearmodulus, and r the material density. A simple solution to this equationis the shear plane wave described by u_(z)=A cos(kx−ct), u_(x)=u_(y)=0.Here c=√{square root over (G/ρ)} is the shear wave velocity and k=w/c isthe wavenumber. The shear modulus can be determined in a material ofknown density through knowledge of the frequency and wavelength of sucha wave. In situations where the shear wave is artificially generated oneof the wave properties may be known a priori. For instance, in“Sonoelastographic imaging of interference patterns for estimation ofthe shear velocity of homogeneous biomaterials,” Physics in Medicine andBiology 49 911-22 (2004), by Wu Z., Taylor L S., Rubens D J., Parker KJ., Wu et al. describe inducing a shear wave of known frequency in themedium and determining the shear modulus by estimating the shearwavelength.

An alternative strategy for determining shear modulus, proposed here, isto create a shear wave of known wavelength and measure its frequency. Ashear wave with particle motion in the z direction generated by a bodyforce f_(z) is governed by the inhomogeneous equation:

$\begin{matrix}{{{\rho \; \frac{\partial^{2}u_{z}}{\partial t^{2}}} - {G\; {\nabla^{2}u_{z}}}} = {f_{z}.}} & (2)\end{matrix}$

Consider the application of an impulsive body force, as from a quickburst of ultrasound, of the form:

f _(z)=δ(t)F(x,y,z)   (3)

in a material initially at rest. The practical utility of the impulse isin approximating a force that occurs on a time scale much shorter thanany other of interest in the problem. In the present problem the forcingfunction pulse can be considered an impulse if its duration is muchshorter than the period of the shear wave it ultimately excites. Forshear waves of on the order of a few kilohertz, this implies a pulseduration of 100 μs or less can be approximated as an impulse.

Treating the force as an impulsive in time, the applied force will existonly at t=0. Since the applied force will be zero for t>0. the shearwave will be governed by the homogenous equation for time t>0, and itsform will be determined by the initial displacement and velocity. Thematerial velocity ∂u/∂t at time 0⁺, the instant after the application ofthe impulsive force, may be determined by substituting f_(z) intoEquation (2) and integrating with respect to time,

$\begin{matrix}{{\int_{- \infty}^{0^{+}}{\left( {{\rho \; \frac{\partial^{2}u_{z}}{\partial t^{2}}} - {G{\nabla^{2}u_{z}}}} \right){t}}} = {\int_{- \infty}^{0^{+}}{{\delta (t)}{F\left( {x,y,z} \right)}{{t}.}}}} & (4)\end{matrix}$

At time t=0⁺, because u_(z)=0 at t=0⁻ and cannot change instantaneously,the G∇²u_(z) term is zero. The velocity at time 0⁺ is then

$\begin{matrix}{{{\rho \; \frac{\partial u_{z}}{\partial t}} = {F\left( {x,y,z} \right)}},} & (5)\end{matrix}$

that is, the initial velocity after the application of the impulsiveforce is proportional to the forcing function.

The propagation of the wave for t>0 may now be calculated using Equation(1) and the initial conditions just determined. For the sake ofsimplicity let us first consider forces and displacements that vary onlywith x and t. In this case Equation (1) simplifies to theone-dimensional equation

$\begin{matrix}{{{\rho \; \frac{\partial^{2}u_{z}}{\partial t^{2}}} - {G\; \frac{\partial^{2}u_{z}}{\partial x^{2}}}} = 0.} & (6)\end{matrix}$

The solution of this equation for t>0 given initial displacement u_(z)=0and initial velocity F(x) is:

${u_{z}\left( {x,t} \right)} = {\frac{1}{2c}{\int_{x - {ct}}^{x + {ct}}{{F(\tau)}{\tau}}}}$

Differentiating with respect to t to obtain the material velocityv_(z)=∂u_(z)/∂t we find that

$\begin{matrix}{{{v_{z}\left( {x,t} \right)} = {\frac{1}{2}\left\{ {{F\left( {x + {ct}} \right)} + {F\left( {x - {ct}} \right)}} \right\}}},\left( {t > 0} \right),} & (7)\end{matrix}$

that is, a wave propagating in the positive and negative x directionwhose form is a replica of the forcing function F(x).

Suppose a forcing function as in Equation (3) with F(x)=1+cos kx isapplied. The material velocity at position x will be given by Equation(7) and simplifies to

v _(z)(x,t)=1+cos(kx)cos(ωt), (t>0),

where

$\begin{matrix}{\omega = {{kc} = {\sqrt{\frac{k^{2}G}{\rho}}.}}} & (8)\end{matrix}$

This velocity signal may be high-pass filtered to remove the DCcomponent.

Measurement of the frequency of oscillation of the velocity functionprovides an estimate of the of the shear modulus G, since the otherquantities in Equation (8) are known. A somewhat more realistic forcingfunction would have the form F(x)=(1+cos kx)ψ(x) where ψ(x) is asmoothly varying envelope function. In this case the velocity signalwould be modulated by the envelope function. A suitable choice ofenvelope function (e.g. Gaussian) will leave the measured frequencyunchanged. The result for an arbitrarily shaped forcing functionF(x,y,z) is not calculated in this work, but simulated using finiteelement methods.

SUMMARY OF THE INVENTION

It is therefore an object of the invention to provide a technique thatprovides a real-time, non-invasive, quantitative measurement of thelocal shear modulus (stiffness) of tissue.

Another object of the invention is to provide a technique that has allthe ARFI advantages but with a quantitative measure of stiffnessincluding when the stiffness is diffuse, i.e., not well defined in anyone area.

To achieve the above and other objects, according to the presentinvention, there is provided a method for determining a shear modulus ofa material, wherein an acoustic radiation force is used impulsively toapply a disturbance of known spatial frequency in a medium of unknownshear modulus and known density. The disturbance gives rise to a shearwave, the velocity and temporal frequency of which are directly relatedto the shear modulus at the initial point of disturbance. Standardultrasonic motion tracking methods may be used to determine thistemporal frequency and thus the local shear modulus. Because thetemporal, rather than spatial, frequency of the shear wave iscalculated, observations of the tissue displacement can be made at asingle point with no need to calculate spatial derivatives ofdisplacement. In addition, calculations of temporal frequency may bemade at several points and averaged to improve the determination, sincethe temporal frequency of the shear wave is fixed by the value of shearmodulus at the point of excitation and does not change with propagationof the shear wave.

The present invention, thus, provides a quantitative measurement of thelocal shear modulus of tissue using a Spatially Modulated AcousticRadiation Force Impulse (SM-ARFI). The SM-ARFI technique in the presentinvention is more advantageous over the known ARFI technique because itincludes all the advantages of an ARFI provides, plus providing aquantitative measurement of the shear modulus of the tissue.

Two papers related to the present invention have been published, whosedisclosures are hereby incorporated by reference in their entiretiesinto the present application:

“Shear modulus estimation by application of spatially modulatedimpulsive acoustic radiation force” by Stephen McAleavey, Manoj Menon,and Jerrod Orszulak, published in Ultrasonic Imaging 29, 87-104 (2007).

“Direct Estimation of Shear Modulus using Spatially Modulated AcousticRadiation Force Impulses”, by Stephen McAleavey and Manoj Menon, 2007IEEE Ultrasonics Symposium, Oct. 28, 2007.

BRIEF DESCRIPTION OF THE DRAWINGS

A preferred embodiment of the invention will be set forth in detail withreference to the drawings, in which:

FIG. 1 shows a flow chart of an exemplary method for determining theshear modulus of a material in accordance with a non-limiting embodimentof the present invention;

FIG. 2 shows a block diagram of an exemplary device that carries out theexemplary method shown in FIG. 1;

FIG. 3( a) illustrates a geometry for Fraunhoffer transducer model;

FIG. 3( b) shows an illustrative model for Gaussian forcing functionsimulation;

FIG. 3( c) shows an illustrative model for linear array forcingfunctions;

FIG. 4 shows SMURF pulse sequence;

FIGS. 5( a)-(c) show graphs of particle velocity vs. time for a uniformmedium with shear modulus G=1 kPa and density r=1000 kg/m³;

FIGS. 5( d)-(f) show graphs of particle velocity vs. time in a phantomwith a G=2 kPa inclusion;

FIG. 6( a)-(c) show graphs of particle velocity vs. lateral position (x)in the direction of the forcing function in a uniform G=1 kPa medium at0 ms, 2 ms, and 10 ms, respectively, after application of the pushingpulse;

FIGS. 6( d)-(f) show graphs of particle velocity vs. position at 0 ms, 2ms and 10 ms, respectively, in a phantom with a G=2 kPa inclusion;

FIGS. 7( a)-(d) show plots of acoustic field intensity generated by 3MHz linear array;

FIG. 8 shows plots of acoustic field intensity generated by 5.33 MHzlinear array;

FIGS. 9( a)-(b) show plots of velocity through time in a uniform G=1 kPaphantom in response to forcing functions generated by 3 MHz focusedlinear array, using the Focal Fraunhofer method at (a) s_(focus)=2 mmand (b) s_(focus)=4 mm;

FIGS. 9( c)-(d) show plots of velocity through time in a uniform G=1 kPaphantom in response to forcing functions generated by 3 MHz focusedlinear array, using the Intersecting Gaussian beams at (c) s_(beam)=2 mmand (d) s_(beam)=4 mm;

FIGS. 10( a)-(h) show plots of velocity vs. position in a uniform G=1kPa phantom in response to forcing functions generated by 3 MHz focusedlinear array, wherein FIGS. 10( a)-(d) are results from using the FocalFraunhofer method with a 2 mm (FIGS. 10( a) and (b)) or 4 mm (FIGS. 10(c)-(d)) standard deviation of the lateral pressure profile. FIGS. 10(e)-(h) are results from using the intersection Gaussian beams with 2mm(FIGS. 10( e) and (f)) or 4 mm (FIGS. 10( g) and (h) standard deviationof the Gaussian describing the beam aperture; and

FIG. 11 shows a plot of velocity vs. time in response to 3 MHz focusedlinear array in the lesion model with G_(lesion)=2 kPa andG_(background)=1 kPa.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENT

A preferred embodiment of the invention will be set forth in detail withreference to the drawings, in which like reference numerals refer tolike elements or steps throughout.

Many non-limiting aspects of the present invention can be betterunderstood with reference to the following figures. The components inthe figures are not necessarily to scale, emphasis instead being placedupon clearly illustrating the principles of the present invention.Further, in describing the non-limiting aspects of the present inventionillustrated in the figures, specific terminology is resorted to for thesake of clarity. Each specific term, however, is meant to include alltechnical equivalents that operate in a similar manner to accomplish asimilar purpose. The present invention, therefore, is not limited to thespecific terms so selected.

FIG. 1 shows an exemplary flow chart of a preferred embodiment of themethod for determining the shear modulus of a material such as human oranimal soft tissue. In step 1, an acoustic radiation force impulse isused to generate disturbance or displacements in a tissue which has aknown density ρ and unknown shear modulus value G. In step 2, thedisturbance at a point in the tissue creates a localized shear wave ofdesired, known wavelength λ and spatial frequency (wave number) k, wherek and λ are related by k=2π/λ. In step 3, the propagation, i.e., shearvelocity c_(s), of the local shear wave is detected by ultrasonic or MRIdevices. In step 4, the temporal frequency of the shear wave iscalculated using the equation:

f _(s) =c _(s)/λ_(s)

It is well known that the shear modulus G of a linear elastic materialsuch as a soft tissue with density ρ is related to the speed of shearwave propagation c_(s) by:

c=√{square root over (G/ρ)}

In terms of the local shear modulus, the temporal frequency of the shearwave is:

$f_{s} = {\sqrt{\frac{G}{\lambda_{s}^{2}\rho}}.}$

From the equation above, the local shear modulus G can be calculatedusing the following equation:

G=ρf_(s) ²λ_(s) ²

In step 5, the local shear modulus G is determined using the knownwavelength λ, known tissue density ρ and determined temporal frequencyf_(s).

Determining the temporal frequency f_(s) is much easier than determiningthe local wavelength of a shear wave. Furthermore, it is well known thatthe temporal frequency of a shear wave remains unchanged even as theshear wave propagates into regions with other values of shear modulus Gdue to the necessity of material continuity during wave propagation withthe spatial frequency k varying to satisfy the wave equation. Thus, thetemporal frequency f_(s) of the shear wave carries information about theshear modulus of the point of excitation. This frequency is unchangedeven if the observation point is located in a region of differing shearmodulus.

Regarding the temporal frequency f_(s), standard ultrasonic motiontracking methods may be used to determine the temporal frequency andthus the local shear modulus. Because the temporal, rather than spatial,frequency of the shear wave is calculated, observations of the tissuedisplacement can be made at a single point with no need to calculatespatial derivatives of displacement.

In FIG. 1, the calculation of temporal frequency in step 4 is made at asingle point in the tissue. However, the calculation of the temporalfrequency may be made at several points and averaged to improve thedetermination because the temporal frequency of the shear wave is fixedby the value of shear modulus at the point of excitation and does notchange with propagation of the shear wave.

FIG. 2 shows an exemplary block diagram of a device 201 that carries outthe method for determining the shear modulus of a soft tissue inaccordance to the preferred embodiment. The device 201 may include anelement 202 such as a transducer to generate acoustic radiation forceimpulses to the soft tissue 203. In response to the disturbance causedby the acoustic radiation force, a shear wave is generated from the softtissue 203. The device 201 further includes a detector 204 that detectsthe temporal frequency of the shear wave from the soft tissue. Thecalculating device 205 calculates the shear modulus using the detectedtemporal frequency, known wave length, and known density of the softtissue as disclosed above. The device further includes an output 206connected to the calculating device to display the value of the shearmodulus of the tissue and other desired values.

A scanner from manufacturer Siemens such as a Siemens Antares may beused as the device 201. The scanner is programmed by the users toprovide custom beam sequences, using a VF7-3 linear array transducer,and to collect the RF echo data. The device 201 can also be anultrasound equipment from companies such as GE and Philips.

Regarding step 1 in FIG. 1 above, the use of acoustic radiation force togenerate displacements in tissue or a medium is known in the art.Acoustic radiation force acts as a body force in an attenuating medium.For plane wave excitation the force is given by F=2aI/c₁, where a is thematerial attenuation coefficient, I the acoustic intensity, and c₁ thelongitudinal speed of sound. This expression has also been found to givequalitatively good agreement with experimental measurements when theintensity is a function of position, as in the beam generated by afocused ultrasound transducer. Therefore by varying the intensity as afunction of position the radiation force may likewise be varied. Linearphased arrays provide great control over the formation of ultrasoundbeams and the acoustic intensity as a function of position. Push beamintensity patterns of a specific lateral spatial frequency may becreated though modulation of the relative phase and/or amplitude ofarray elements. Several methods are described in the literature forcreating such a pattern. We describe four of these methods here:

The first is the Focal Fraunhofer approach. It is well known that in thefar field of an unfocused transducer or in the focal zone of a focusedtransducer the lateral pressure variation is the Fourier transform ofthe apodization function of the aperture,

p(x,y,z)=P ₀λ² z ²∫∫_(S) A(λzx ₀ , λzy ₀)e ^(−f2πxx) ⁰ e ^(−j2πxy) ⁰ dx₀ dy ₀

where A(x, y) is the complex apodization including the phase shift forfocusing and

$P_{0} = \frac{j\; \rho \; c_{l}v_{0}^{{j\; 2\pi \; f_{l}t} - {k_{l}z}}^{{- j}\; {{k_{l}{({x^{2} + y^{2}})}}/2}z}}{\lambda_{l}z}$

where v₀ is the peak velocity amplitude of the transducer, and k₁, c₁,and l₁ are the wavenumber, sound speed, and wavelength of thelongitudinal (i.e. ultrasound) wave. The aperture geometry isillustrated in FIG. 3 a. For separable A(x, y)=A_(x)(x) A_(y)(y) thepressure in the y=0 plane is simply

p(x,z)=P ₀ P _(y)λ₁ z∫A(λ₁ zx ₀)e ^(−j2πxx) ⁰ dx ₀   (9)

where P_(y)=∫A_(y)(y₀)dy₀. Therefore a sinusoidal oscillation ofpressure may be induced by an apodization consisting of a pair ofimpulses. This method was employed by Anderson to generate a laterallyvarying pressure for estimation of the lateral component of blood flowvelocity. In order to radiate any significant energy, an apodizationconsisting of a pair of well-localized windows is used. For instance, anapodization consisting of a sum of Gaussians displaced from the beamaxis of the form

A(x ₀)=e ^(−(x) ⁰ ^(−d)) ² ^(/2s) ² +e ^(−(x) ⁰ ^(+d)) ² ^(/2s) ²   (10)

gives a pressure at the focus z_(f) of the form

p(x)=P ₀ P _(y)2s√{square root over (2π)}e ^(−(kxx)) ² ^(/2z) ^(f) ²cos(kxd/z _(f))   (11)

and an intensity of

I(x)=I ₀4πs ² e ^(−(ksx/z) ^(f) ⁾ ² (1+cos(2kxd/z _(f)))   (12)

The limitation of this technique can be seen in Equation (12). As s isincreased, allowing more energy to be put into the beam, the focalregion narrows proportionately. This relationship makes it difficult tocreate a high pressure in the focal zone while also generating a lateralpressure distribution of well-defined spatial frequency. For theintensity signal

${I(x)} = {I_{0}4\pi \; \frac{z^{2}}{2k^{2}\sigma^{2\;}}{^{{{- x^{2}}/2}\sigma^{2}}\left( {1 + {\cos \left( {2\pi \; {x/\lambda_{s}}} \right)}} \right)}}$

the apodization parameters are

$d = {{\frac{z_{f}\lambda_{0}}{2\lambda_{s}}\mspace{14mu} {and}\mspace{14mu} s} = {\frac{\sqrt{2}\lambda \; z_{f}}{4\pi \; \sigma}.}}$

A second method, the intersecting plane wave approach, is to useintersecting unfocused beams of identical frequency traveling at anoblique angle to each other. The simplest expression of this idea is apair of plane waves propagating at angles of ±θto the z axis,

p(x,t)=p ₀ cos(2πf ₁ t+k ₁(z cos θ−x sin θ))+p ₀ cos(2πf ₁ t+k ₁(z cosθ+x sin θ))

which is equal to

p(x,t)=2p ₀ cos(2πf ₁ t+z cos θ)cos(k ₁ x sin θ).

Such plane waves set up a lateral pressure variation of wavelengthl_(s)=1₁/sin q. The intensity is given by

I(x)=4p ₀ ²(1+cos(2k ₁ x sin θ))/2   (13)

A finite aperture cannot be used to generate plane waves, but a goodapproximation can be achieved over a limited area by intersecting twounfocused Gaussian beams. Such beams maintain their lateral profileuntil the near/far transition, when they begin to spread. The advantageof this method compared to the Focal Fraunhofer method is that, whilethere is no focal gain, a greater number of elements may be used to formthe beams while still creating a wide region of excitation.

A third method, the sequential excitation, is to generate the desiredintensity distribution sequentially, rather than simultaneously, byfiring a number of beams with intensity and foci located so that the sumintensity over the entire sequence is the desired intensity. This methodis feasible because the great difference between longitudinal ultrasoundwave speed (1540 m/s) and shear wave speed (1-10 m/s) and highultrasound pulse repetition frequencies that may be used allow an entirebeam sequence to be transmitted before the shear wave has propagatedappreciably. Indeed, since one need not wait for echoes from the pushbeam, the PRF may be as high as the reciprocal of the push time at anyone location regardless of the desired depth. The advantage of thismethod is that much greater intensities can be achieved, because a largenumber of transducer elements may be used to generate each beam in thesequence, unlike the Focal Fraunhofer method, while at the same timebeam focusing is implemented to achieve high local intensities, unlikethe intersecting plane waves method.

In a fourth method, we apply the constrained least-squares method ofWalker and Guenther to determine apodization and phase vectors in termsof λ₁, λ_(s), and ν₀. Simultaneously focusing several parallel beamsusing interlaced apertures, we also consider the generation of multiplelateral foci rather than a pair of intersecting beams. The focusingmethods considered will include simple geometric focusing, Gaussian beamphasing/apodization, and interlaced non-diffracting (X-wave) apertures.The results will be evaluated both in terms of maximizing the peakintensity generated by each method in the region of interest andminimizing the intensity outside it.

Of the four methods, the Focal Fraunhofer and Gaussian beam methods aredemonstrated in this paper. The result of any of these methods is alateral modulation of beam intensity of the form I(x)=(1+cos k_(s)x)ψ(x)where ψ(x) is a smoothly varying non-negative envelope function. It isimportant to remember that the spatial frequency of the lateralvariation in intensity, k_(s), is twice the spatial frequency of thepressure variation, due to the squaring of pressure to calculateintensity. The shear wave generated will be of spatial frequency k_(s).

Simulation will now be described. Two-dimensional finite element modelswere employed to model the response of tissue to spatially varyingacoustic radiation force and the propagation of shear waves from theexcitation location. Comsol Multiphysics (Comsol, Inc.) was used tocreate and analyze the 2D plane-strain model. A mesh pitch of 0.25 mmwas used for all models to provide adequate spatial sampling. Theadequacy of this pitch was verified by mesh refinement tests. Tissue wasmodeled as a linear elastic solid with a density of 1000 kg/m³.Radiation force was modeled as a body force in the direction of beampropagation. Two models were created.

The first model, Gaussian Force Model, depicted in FIG. 3 b, was auniform elastic solid of dimensions 6 cm×6 cm (the equivalent 3D solidhaving identical properties and forces that do not vary perpendicular tothe cross-section shown). The solid was modeled as having a density of1000 kg/m³, a Poisson's ratio of 0.49, and shear modulus G of 1 kPa.This model was evaluated for Young's modulus values of 3 kPa and 6 kPa(G=1 and 2 kPa). A Raleigh damping model with a=1 and b=0 was used.Lagrange quadratic elements were used. A fixed boundary condition wasapplied to prevent bulk motion of the model. In this model the radiationforce was modeled as a Gaussian-weighted cosine applied for 50 μs anddescribed by

F _(z)(x,z)=F ₀(1+cos 2000πx)e ^(−(x) ² _(+y) ² ^()/(2·0.001) ² ⁾   (14)

where F₀ is the peak force. This forcing pattern had a spatial frequencyin the x direction of 1 cycle/mm, corresponding to a shear wavefrequency f_(s) of 1 kHz in 1 kPa shear modulus material.

The second model, the Linear Array Model, used forcing functionsgenerated by linear arrays using techniques described in the Descriptionof Related Art section above. In this model the transducer was centeredat the origin with the z axis normal to the transducer face. The modelwas either uniform with a shear modulus of 1 kPa, or included a 0.8 cmdiameter 2 kPa lesion centered at a depth of 4 cm. The Poisson's Ratiofor this model was 0.4995, and the Raleigh damping coefficients were α=1and β=10⁻⁵. This slightly larger value of b improves the numericalconvergence and runtime of the finite element program and does notsignificantly alter the simulation results. A more realistic model ofthe applied radiation force was used for this model. The Field IIultrasound simulation program was used to calculate beam patterns fortwo linear arrays. The arrays had center frequencies of 3 MHz and 5.33MHz and were simulated using the parameters given in Table 1 below. Thesimulation geometry is depicted in FIG. 3 c.

TABLE 1 Transducer Parameters for Field II Simulations Transducer ATransducer B Frequency 3.0 MHz 5.33 MHz Element Pitch 0.4 mm 0.2 mmElement Height 15 mm 7.5 mm Elevation Focus 4 cm 3.75 cm

Both the Focal Fraunhofer and intersecting Gaussian beam methods wereexamined. In all cases a lateral intensity modulation of 1 cycle/mm wasgenerated. In the Focal Fraunhofer method simulation the element phasingwas selected to form a focus at z=4 cm with transducer A and at z=2 cmwith transducer B. The apodization was determined using Equation (12) togenerate a cosine-modulated Gaussian intensity profile at the focus. Forarray A apodizations were calculated for obtaining standard deviationsof the Gaussian lateral profile of s_(x)=2 mm and s_(x)=4 mm, while forarray B a focal width of s_(x)=1 mm was selected. The intersecting beammethod was implemented by generating pairs of Gaussian apertures withstandard deviations of s=2 mm and s=4 mm in the case of transducer A ands=1 mm for transducer B. The phase across each aperture was variedlinearly to steer the beams by ±q and generate the desired lateralpressure modulation spatial frequency, as given by Equation (13). Thedistance between the aperture pair centers was selected so the beam axesintersected at z=4 cm for transducer A and z=2 cm for transducer B.

The use of Field II allows comparisons between the intensities achievedby the different beamforming methods. This is especially useful fordetermining which method will generate the largest intensities (andhence displacements) given a limited maximum drive voltage for eachtransducer element. A tissue attenuation value of 0.7 dB/cm/MHz was usedand non-linear propagation effects were not considered. Intensity valuesfor a 50 μs pushing pulse were calculated on a 10 points/mm grid in thescan (y=0) plane. The intensities at each point were converted to a bodyforce using the F(x, y)=2aI(x, y)/c relationship with the direction ofthe force away from the transducer.

In vitro experiments were performed with a Siemens Antares ultrasoundscanner

(Siemens Medical Solutions USA, Ultrasound Group). A VF7-3 linear arraytransducer, which roughly matches the characteristics of Transducer B inTable 1 in FIG. 12 was used. A custom beam sequence, illustrated in FIG.4, was created to generate a modulated radiation force and track theresulting shear waves. The beam sequence is similar to the shear wavetracking sequence described by Nightingale K., McAleavey S., and TraheyG. in “Shear-wave generation using acoustic radiation force: In vivo andex vivo results,” Ultrasound in Medicine and Biology 29, 1715-23,(2003). The pushing beam is fired at the same location for allensembles, while the tracking sequence is swept across the region ofinterest. By repeatedly firing the pushing pulse at a single locationand tracking at multiple locations a composite sequence is developedwhich tracks the progress of the shear wave through the entire region ofinterest. This sequence allowed the propagation of the induced shearwave to be visualized. However, such visualization is not necessary fordetermining frequency of the induced shear wave. The shear wavefrequency and shear modulus may be determined from a single excitationand tracking of a single point.

Tracking and pushing beams were both focused at a depth of 2 cm. Thetracking beams employed a raised-cosine apodization. The pushing beamwas generated using a pair of rectangularly-apodized apertures withlinear phasing across each to steer the beams to intersection andgenerate a 1 cycle/mm lateral intensity variation. The pushing pulse was200 cycles long with a center frequency of 6.67 MHz. A total of 64elements were used, 32 each for the left and right apertures.

Tissue mimicking phantoms, similar to those described by Hall T. J.,Bilgen M., Insana M. F., and Krouskop T. A. in “Phantom materials forelastography,” IEEE Transactions on Ultrasonics, Ferroelectronics andFrequency Control, 44, 1355-65, (1997), were used for in vitro testingof the method. Two phantoms were constructed, one “soft” and one“stiff”. The soft phantom had a gelatin concentration of 33 g/literwater, while the stiff phantom had a concentration of 99 g/liter water.Graphite powder at a concentration of 160 grams/liter was added to actas a scattering agent, n-propanol was added at a concentration of 50ml/liter to adjust the longitudinal sound speed of the phantom, andgluteraldehyde at 10 mL/liter to cross-link the gelatin. The phantomswere nearly cylindrical, with a top diameter of 6.5 cm, bottom diameterof 8 cm, and height of 7 cm.

The shear modulus of each phantom was determined using unconfined cyclicloading to estimate Young's modulus. The shear modulus can be found fromthe Young's Modulus E using the relationship G=E/2(1+n), where n isPoisson's ratio, assumed here to be 0.5 (incompressible). The entirephantom was placed in between two oil-coated platens of a MTS SintechMaterial Testing Workstation (MTS System Corporation). The phantomsunderwent compression up to a 5% strain at a rate of 2 mm/s while theforce was measured using a 5N load cell. Young's modulus was determinedfrom the slope of the stress-strain curve in the 3-5% strain region toavoid measurement error due to contact non-uniformity between thephantom and platens.

For the results of the Gaussian Model simulation, the velocity versustime for the Gaussian forcing function of Equation (14) in a uniform G=1kPa function at the points (x,z)=(0,0, (3,0), and (7,0)mm is shown inFIG. 5 a-c. In all cases the time between wave peaks is 1 ms, asexpected from Equation (8). Slight disturbances due to the reflection ofa longitudinal wave are visible in the plots. These are discernible dueto the relatively low Poisson's ratio (0.49) used in the model. Valuesof Poisson's ratio closer to 0.5 would diminish the amplitude of thelongitudinal wave.

Graphs of velocity versus time for the same Gaussian forcing functionacting at the center of a G=2 kPa region 1 cm² enclosed in a G=1 kPabackground (see FIGS. 3 a-c) are presented in FIG. 5 d-f at the samespatial locations as plots a-c. We note that the doubling of shearmodulus has resulted in a √{square root over (2)} increase in thetemporal frequency of the wave, estimated from the peak-peak interval,as expected from Equation (8). The peak of the wave arrives at the 3 mmpoint at 2.1 ms in keeping with the 1.4 m/s velocity within the lesion.Notice also that the temporal frequency of the wave does not change asit propagates from the G=2 kPa region to the G=1 kPa region, asexpected. The peak arrives at the 7 mm point at 5.6 ms, as expected forpropagation through 5 mm at 1.4 m/s and 2 mm at 1 m/s.

FIG. 5 reveals that the velocity wave generated by the Gaussian forcingfunction contains negative velocity values. The one dimensional model ofEquation (7) suggests that, since the forcing function is non-negative,the velocity wave should likewise be zero or greater. The discrepancy isdue to the the fact that the forcing function used varies both in x andz, and the Laplacian in Equation (1) will contain a non-zero ∂²u_(z)/∂z²term. In the one-dimensional case the forcing function and displacementdo not vary with z and the ∂²u_(z)/∂z² is zero. In both cases theessential aspect of the forcing function, the imposed wavelength, ismaintained.

Velocity as a function of x in the uniform model is plotted in FIGS. 6a-c at t=0, 2 and 10 ms. The wave propagating through a uniform G=1 kPamaterial maintains its 1 mm wavelength. In FIGS. 6 d-f the contractionof wavelength by a factor of 1/√{square root over (2)} as the wavepropagates from the G=2 kPa (x<5 mm) to the G=1 kPa (x>5 mm) is visible.Also visible in FIG. 6 f at about 4 mm is the a reflected wave from theleft-hand boundary of the lesion. That initial velocity matches theforcing function can be seen clearly in FIGS. 4 a and d. The initialvelocity is also insensitive to the modulus of the lesion. This followsfrom Equation (4), where the shear modulus multiplies only the Laplacianof displacement, which is zero at time zero.

For the results of the linear array model simulation, the intensitypatterns generated by the 3 MHz and 5.33 MHz linear arrays as calculatedby Field II are plotted on a linear grayscale in FIGS. 7 a-d and FIG. 8.For clarity each plot is normalized to the maximum calculated intensityin the field of view. The fields generated by the focal Fraunhofermethod (FIGS. 7 a and b, FIG. 8 left) show a pronounced broadening withincreasing distance from the transducer. This leads to a depth-dependentexcitation wavelength, which would cause a bias in the shear modulusestimates were it not accounted for. The intersecting Gaussian beammethod shows much smaller broadening, since the region of interest(around z=4 cm for the 3 MHz array and z=2 cm for the 5.33 MHz array) isin or close to the near field of the individual Gaussian apertures.

Another significant difference between the two forcing functiongeneration methods is in the relative near-field intensities. Comparedto the intersecting Gaussian beams methods, the near-field intensitiesgenerated by the Focal Fraunhofer method are greater in proportion tothe peak intensity generated at the focus. This result shows up thelimitation of the Focal Fraunhofer method, namely that a small number ofelements may be used to generate the desired lateral intensitymodulation, and that as a consequence a greater proportion of thepush-beam energy will be lost to the near field.

The response of the uniform G=1 kPa phantom to the forcing functions ofFIG. 7 is shown in FIGS. 9 a-d as a function of time and in FIG. 10 as afunction of lateral position. In all cases the greater force exerted bythe intersecting Gaussian beams method is evident, increasingly so asthe width of the region of excitation is increased. The initialvelocities are observed to be a good match to the lateral distributionof the beam intensities. The response through time at two points((0.4)cm and (0.5,4)cm) near the “focus” of the pushing beams is shownin FIGS. 9 a-d. In all cases the expected 1 kHz frequency is observed.The 1 kHz signal rides on top of a low-frequency component due to theaxial and lateral variation in the envelope of the forcing function,which may be removed by high-pass filtering. The evolution of the shapeof the velocity wave as it propagates from its initial non-zero form isillustrated in FIG. 10, where velocity as a function of position isplotted at 0 and 5 ms after the application of the forcing function.

The velocity through time in the 0.8 cm lesion model is shown in FIG.11. The frequency of the velocity wave at the center of the lesion andat a point 2 mm outside are seen to vary by the expected factor of√{square root over (2)}. This reassuringly illustrates that, though theforced region is significantly larger than the lesion, it is possible todistinguish the lesion from the background by measuring the localvibration frequency. Results obtained with the other pushing beams weresimilar.

For the in vitro results, the two gelatin phantoms had Young's modulusvalues of 16.9±0.5 kPa and 4.6±0.2 kPa, as estimated by the MTSmeasurement system. These correspond to shear moduli of 5.6 kPa and 1.5kPa, determined by the relationship G=E/2(1+n), where we have assumedthe phantom is incompressible and Poisson's ratio n=0.5. Theexperimentally measured velocity though time in response to the 1cycle/mm pushing beam in the two gelatin phantoms. The value of shearmodulus for each phantom is determined using Equation (8), with k=2p/l=2000 pm⁻¹ and density r=1000 kg/m³. The period of oscillation of thesoft phantom is estimated to be 0.83 ms (based on peak-to-peak time),corresponding to a frequency of 1.2 kHz and a shear modulus of 1.4 kPa.The stiff phantom has an estimated modulus of 5.8 kPa. Both values showgood agreement with the MTS measurements.

The success of the method in the present invention depends on theability to generate a nearly impulsive forcing function with awell-defined spatial frequency. The approximation of the forcingfunction is well justified for pulses on the order of 10-100 μs, sincethe induced wave will only propagate a small fraction of a wavelengthduring the excitation period. For instance, a l=1 mm wave in a G=2 kPamedium will propagate only 71 μm or 1/14 of a wavelength over theduration of a 50 μs pushing pulse. Two methods have been demonstratedhere based on either intersecting Gaussian beams or the Fraunhoferapproximation. While either method can generate the desired spatialfrequency, for a given peak pressure from any given transducer element,the Gaussian beam method is able to deliver a higher intensity to theregion of interest for the model considered here, due to the greaterutilization of elements in the aperture. In “Supersonic shear imaging: Anew technique for soft tissue elasticity mapping” IEEE Transactions onUltrasonics, Ferroelectronics and Frequency Control 51, 396-409, (2004),Bercoff J., Tanter M., and Fink M. have demonstrated the synthesis ofspatially extended shear wave sources using radiation force from asequence of pushing beams. It may be feasible to create spatiallymodulated beams of well defined spatial frequency using this method aswell. The pulse-repetition frequency used for sequentially “painting”the spatially modulated pushing beam may be very high, since there is norequirement to wait for echoes to return from any depth. Therefore, theeffective pulse repetition period could be on the order of the durationof an individual pushing beam, 50 μs, creating a nearly impulsive forcefrom a sequence of pushes.

The analysis presented here assumes that propagation of shear wavesthrough the medium to be interrogated is purely elastic and lossless.Tissue is known to be viscoelastic. Among the consequences ofviscoelasticity are attenuation and frequency dependent propagationvelocity (dispersion) of shear waves. In comparison withsonoelastography, which typically uses vibration frequencies in the100-500 Hz range, the frequencies of the shear waves excited in thesimulations and experiments here are higher, on the order of one to afew kilohertz. Higher and lower frequencies can be generated easily byvarying the spatial frequency of the pushing beam. Higher frequencyshear waves will be attenuated more rapidly, but because they aregenerated by acoustic radiation force at the point of interest, theyneed propagate no more than a few wavelengths to enable a measurement.Therefore, tissue may be examined over a range of frequencies where theuse of externally generated shear waves would be precluded byattenuation. The application of spatially modulated radiation force asdescribed gives rise to a shear wave with a well defined frequency, andpoints to a potential to measure the frequency dependence of shear wavespeed. In a simple elastic medium, the applied wavelength and measuredfrequency are linearly related. In a viscoelastic medium where shearvelocity is a function of frequency, we would expect a deviation fromthis linear behavior. By applying a sequence of spatially modulatedpushing pulses over a range of spatial frequencies and measuring theinduced frequency the wave speed at the observed frequency may beestimated.

In conclusion, a method for the calculation of shear modulus usingspatially modulated impulsive acoustic radiation force has beenpresented. In both the simulations and in vitro experiment we find goodagreement between the modulus estimated by spatially modulated radiationforce and ground truth or MTS measurements. The advantage of theproposed technique compared with other methods is the relative ease withwhich vibration frequency may be estimated within tissue, compared tothe difficult problem of estimating shear wavelength in tissue. Theproposed method requires at most high-pass filtering the ultrasonicallymeasured displacement to remove low-frequency motion. In contrast,inversion techniques often require the estimation of second derivativesof displacement with respect to space and time, resulting in a verynoise-sensitive estimate.

The foregoing description and drawings should be considered asillustrative only of the principles of the invention. The invention maybe configured in a variety of shapes and sizes and is not intended to belimited by the preferred embodiment. Numerous applications of theinvention will readily occur to those skilled in the art. Therefore, itis not desired to limit the invention to the specific examples disclosedor the exact construction and operation shown and described. Rather, allsuitable modifications and equivalents may be resorted to, fallingwithin the scope of the invention.

We claim:
 1. A method for determining a shear modulus of a material,comprising the steps of: applying a disturbance having a spatialfrequency to the material at an initial point, causing a shear wave tobe produced in the material; detecting the shear wave; determining atemporal frequency of the shear wave; and determining the shear modulusof the material based on the determined temporal frequency.
 2. Themethod as recited in claim 1, wherein the material comprises human oranimal soft tissues.
 3. The method as recited in claim 1, furthercomprising determining temporal frequencies at several points of thematerial and averaging to determine an average temporal frequency forthe step of determining the shear modulus of the material.
 4. The methodas recited in claim 1, wherein the step of determining the temporalfrequency comprises determining the temporal frequency throughultrasonic motion tracking.
 5. The method as recited in claim 1, whereinthe material has a known density, and the step of determining the shearmodulus of the material is also based on the known density and awavelength of the spatial frequency.
 6. The method as recited in claim1, wherein the step of applying the disturbance having the spatialfrequency comprises using an acoustic radiation force impulse togenerate the disturbance.
 7. An apparatus for determining a shearmodulus value of a material having a known density, comprising: atransducer for generating at least one acoustic radiation impulse towardthe material, the at least one acoustic radiation impulse having knownspatial frequency and wavelength to produce a shear wave propagating inthe material; a detecting device for detecting the shear wave anddetecting a temporal frequency of the shear wave; a calculating devicefor calculating the shear modulus of the material using the detectedtemporal frequency, known density of the material, and known wavelengthof the spatial frequency; and an output coupled to the calculatingdevice for indicating the shear modulus value of the material.
 8. Anapparatus as disclosed in claim 7, wherein the calculating devicedetermines temporal frequencies at several points of the material andaverages to determine an average temporal frequency.
 9. An apparatus asdisclosed in claim 7, wherein the detecting device detects the temporalfrequency through ultrasonic motion tracking.